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ABSTRACT 

Adaptive  optics  systems  are  designed  primarily  for  astronomical  and  surveillance 
applications  to  compensate  for  the  effect  of  random,  time-varying  aberrations  introduced  by  the 
earth's  atmosphere.  Compensation  for  these  aberrations  predicates  some  means  of  sensing  the 
atmospheric  distortion  with  sufficient  signal-noise  ratio  to  apply  a  meaningful  correction.  In  most 
scenarios  of  interest  (e.g.  remote  stellar  objects,  earth-orbiting  satellites),  the  light  level  available 
for  this  process  is  low,  placing  limits  on  the  degree  of  correction  that  is  possible.  The  use  of 
artificial  guide  stars  (laser  beacons)  is  one  approach  to  improving  the  signal-noise  available  for 
wavefront  sensing.  However,  the  technical  effort  and  cost  associated  with  this  solution  is  very  high 
and  still  technical  problems  remain  to  be  solved  (e.g.  the  tilt  problem).  Certain  surveillance 
applications  may  also  disqualify  solutions  which  employ  target  illumination.  There  is  thus 
motivation  for  the  use  of  systems  using  simpler  technology  but  employing  optimally  configured 
sensing  to  maximise  performance. 

The  overall  study  presented  is  concerned  with  the  development  of  computer  simulation 
tools  enabling  the  study  of  optimum  imaging  scenarios  for  low  light  level  adaptive  optics  and  of 
assessing  their  implications  for  optical  interferometers  in  which  the  constituent  telescopes  employ 
adaptive  optics.  The  first  part  of  the  report  outlines  the  overall  scope  and  scientific  basis  of  an 
adaptive  optics  simulation  program  developed  during  the  study.  This  includes  a  description  of  a 
graphical-user-interface  (GUI)  which  has  been  developed  to  assist  the  easier  use  of  the  program  and 
some  explanatory  material  of  the  current  version  of  the  simulation  program.  The  second  part  of  the 
report  details  results  of  some  low  light  level  studies  on  a  lm  telescope  employing  adaptive  optics 
without  laser-guide  star  technology.  The  implications  for  low-light  level  AO  applications  are 
briefly  discussed.  Part  3  is  given  in  the  form  of  a  journal  article  written  in  the  context  of  this  study 
and  submitted  to  J.O.S.A.  letters  entitled  -  “Variational  solution  for  minimum  error  norm 
wavefront  projectors”.  The  paper  describes  the  derivation  of  an  optimal  method  for  estimating 
atmospherically  distorted  wavefronts  by  direct  integration  of  the  wavefront  slope  measurements 
which  is  also  guaranteed  to  have  minimum  noise  propagation. 
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1  PHYSICAL  BASIS  FOR  ADAPTIVE  OPTICS  COMPUTER  MODEL 

1.1  Software  modules 

The  adaptive  optics  computer  model  uses  6  basic  software  modules  which,  interfaced  in  a 
straightforward  manner,  simulate  the  behaviour  of  an  adaptive  optics  system.  The  modules 
correspond  to  - 

•  The  generation  of  a  sequence  of  time-evolving  atmospheric  wavefronts. 

•  Image  formation  by  a  Shack-Hartmann  sensor  coupled  to  a  CCD  camera  and  estimation  of 
atmospheric  wavefront  slopes  by  sub-aperture  centroiding. 

•  Modal  estimation  of  the  atmospheric  wavefronts. 

•  Modelling  of  the  spatio-temporal  behaviour  of  the  adaptive  mirror  and  thus  the  drive  voltage 
transfer  matrix. 

•  Implementation  of  a  closed-loop  control  algorithm,  capture  of  data  in  the  image  plane  of  the 
telscope  and  generation  of  appropriate  performance  indicators  such  as  time-averaged  mean- 
square  error  or  Strehl  ratio. 

•  A  graphical  user-interface  to  allow  easy  use  of  the  program. 

The  scientific  basis  of  the  computer  model  is  generally  accurate  and  realistic  but  certain 
approximations  have  been  made.  Accordingly,  we  first  outline  the  basic  principles  underpinning 
our  computational  implementation  of  these  modules. 

1.2  Generation  of  time-evolving  atmospheric  wavefronts 

The  time-evolving  wavefront  generator  is  based  on  the  cumulative  effect  of  4 
atmospheric  layers,  each  modelled  by  a  random  phase  screen.  The  phase  power  spectrum  is 
modelled  by  the  Von-Karman  spectrum. 

The  basic  approach  to  simulation  is  based  on  random  sampling  for  the  phase  as  a  function  of 
space  and  time  such  that  it  will  statistically  obey  the  required  space-time  phase  correlation 
function  [1],  defined  as  - 

rv  =  (<p{Xi,t,)<p(xj,tm))  t1  1] 

Where  ^(x,,?*)  represents  the  phase  evaluated  at  point  x,  and  time  tk. 

We  represent  the  simulated  wavefront  by  an  array  of  spatial  phase  values  evaluated  at  N  co¬ 
ordinates  [x,  ].=1  v  over  the  pupil.  Further,  let  us  assume  that  we  want  to  generate  M  such 
successive  phase  screens.  This  sequence  of  phase  screens  can  be  represented  by  a  single  vector  a 
of  length  NxM: 
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We  can  then  write  the  space  time  covariance  matrix  of  the  phase  screens  as  the  covariance  of 
that  vector  (a  matrix  of  size  (NxM)2): 

r,=r.=(a„a;>  1*31 


The  correlation  function  of  the  phase  is  related  to  the  structure  function  of  the  phase  and  can  be 
evaluated  for  a  given  model  of  the  refractive  index  power  spectrum  such  as  the  Von  Karman  or 
Kolmogorov.  However,  the  Kolmogorov  spectrum  presents  a  singularity  at  the  spatial  frequency 
zero.  It  can  be  shown  that  the  covariance  matrix,  calculated  from  this  spectrum,  will  also  contain 
such  a  singularity.  If  we  considered  a  piston  removed  phase,  i.e.  a  phase  function  with  a  mean  of 
zero,  the  covariance  matrix  does  not  contain  any  singularity  but  still  results  in  an  extremely 
complex  expression  for  the  covariance  matrix  [1].  This  singularity  does  not  occur  if  we  assume 
the  Von-Karman  spectrum.  The  covariance  matrix  corresponding  to  the  Von-Karman  spectrum 
is  real,  symmetric  and  positive  definite.  A  Cholesky  factorisation  of  that  matrix  into  two  square 
matrices  is  therefore  possible: 

Ta  =  RRT  [14] 

The  method  proceeds  by  generating  a  vector  of  length  NxM  of  independent  gaussian  random 
variables  b  and  a  product  vector  c  =  Rb  from  which  it  follows  that  - 

co  v(c)  =  (ccT  ^ 

=  {RbbTRT) 

=  R{\}bT^Rr  (As  R  is  not  a  random  quantity) 

[1.5] 

The  components  of  the  Gaussian  vector  b  are  statistically  independent  which  means  that  there  is 
no  correlation  between  them.  Further,  we  may  choose  them  to  have  unit  variance  so  that  the 
covariance  reduces  to- 

cov(c)  =  RRT  =  [16] 

as  required. 


We  have  shown  above  that,  by  random  generation  of  a  vector  of  unit  variance  gaussian 
variables,  we  can  obtain  a  set  of  random,  temporally  correlated  phase  screens  obeying  the 
correct  spatio-temporal  statistics.  This  method,  however,  does  require  significant  computational 
resources.  Consider,  for  example,  if  we  wish  to  generate  five  successive,  correlated  phase 
screens  of  size  128x128  pixels.  This  implies  the  need  to  compute  a  covariance  matrix  of  size: 
(128xl28x5)2=  6710886400  «  6.7  *109  elements.  Cholesky  factorisation  of  such  a  matrix,  even 
by  today’s  standards,  requires  a  very  powerful  computer.  To  partly  circumvent  this  problem,  we 
first  simultaneously  generated  only  three  consecutive  screens.  Then,  by  varying  the  time  step 
between  each  screen  in  this  set,  we  may  subsequently  reconstruct  the  required  sequence  of  time- 
evolving  phase  screens  by  a  temporal  interpolation  procedure.  For  example,  if  we  want  a 
sequence  of  phase  screens  over  a  time  interval  of  length  T=2s,  with  a  time  separation  between 
each  screen  of  r  =  2ms  we  first  generate  a  stack  of  three  wavefronts  separated  by  a  time  step 

—  -  r  Then  we  successively  decrease  the  time  step  until  we  produce  all  the  necessary  screens. 
2 

To  further  conserve  computational  resources,  we  may  decrease  the  spatial  size  of  each  screen, 
evaluating  the  phase  values  at  only  a  small  number  of  selected  spatial  co-ordinates.  Provided  we 


T 

use  a  sampling  interval  smaller  than  y,  then  we  can  spatially  extrapolate  the  randomly 

generated  phase  screens  between  the  sampled  points.  Details  of  the  technique  can  be  found  in 
reference  [1], 

1.3  Shack-Hartmann  image  formation  and  centroiding 

The  Shack-Hartmann  sensor  has  been  widely  studied  (see  reference  [2]  for  an  exhaustive  list  of 
references)  and  is  the  most  common  form  of  wavefront  sensor.  This  is  largely  due  to  it  s  conceptual 
simplicity  and  competitive  performance.  In  essence,  a  Hartmann  sensor  provides  a  number  of 
wavefront  slope  estimates,  each  averaged  over  a  sub-region  of  the  pupil,  by  calculation  of  the 
barycentre  of  each  sub-aperture  image.  Our  implementation  of  Shack-Hartmann  image  formation  is 
a  high-fidelity  model  and  includes  realistic  effects  of  low  light  level  imaging  including  both 
detector  and  photon  noise.  The  details  of  the  method  are  outlined  in  Ash  et.  al.  [3],  The  centroiding 
method  employed  in  our  model  is  based  on  a  recursive  centre-of-gravity  technique.  It  has  been 
shown  that  matched  filtering  techniques  [4]  may  offer  improved  slope  estimates  at  low  light  levels 
but  the  current  version  of  the  simulation  program  does  incorporate  such  an  estimator. 

1.4  Modal  estimation  of  atmospheric  wavefronts 

This  subject  has  been  widely  studied.  We  briefly  summarise  our  approach.  The  basic 
method  is  to  expand  the  random  wavefront  in  terms  of  a  set  of  chosen  basis  functions  ■?*(*)- 

<p{x)  =  YjakPk^) 

Since  the  Hartmann  sensor  provides  a  finite  set  of  local  phase  gradient  measurements,  we  take  the 
gradient  of  this  equation  and  evaluate  at  the  coordinates  (x,.)at  which  the  gradient  is  known  - 

V^(x,)  =  X^VPi(xj)  (i  =  1,2.. .M)  [1.8] 

This  system  may  be  written  in  matrix  form  as  - 

Ax  +  e  =  m  [19] 

We  seek  a  linear  estimator  for  the  unknown  modal  coefficients  given  by  a  linear  combination  of 
the  slope  data  - 

i  =  Lm  =  L(Ax  +  e)  [110] 

which  will  minimise  the  error  covariance  matrix  given  by  - 
E  =  ^(x-x)(x-x)T^  [111] 

Where  the  matrix  L  is  to  be  determined.  Inserting  [1.10]  into  [1.11]  and  using  the  definition  of 
Cx  =  (xxT  } ,  we  obtain  - 

E  =  P-LA]Cx[I-LA]t+LVLt  [1-12] 

Taking  first  variations  in  E  with  respect  to  L  and  setting  to  zero  yields  - 
<SE  =  <5L(ACX[I  -  LA]t  +  VLt  )  +  ([I  -  LA]CXAT  +  LV)<SLT  =  0  [1.13] 

The  optimal  matrix  is  thus  given  by  - 


Lmap  -  CxAt(ACxAt  +  v) 


[1.14] 


This  may  also  be  written  as  - 

LOFr=(ATV-lA  +  C;1)_1ATV-1  [1.15] 


Substituting  this  optimal  (MAP)  estimator  into  the  definition  of  the  error  covariance  matrix  (eq. 
1.10)  gives  the  expression  - 


EOFr=(ATV-1A  +  Cl-1) 


[1.16] 


•  If  we  assume  that  we  have  no  prior  knowledge  of  the  modal  coefficients  (i.e.  Cx  — »  oo ),  this 

optimal  estimator  reduces  to  the  well-known  BLUE  estimator  - 

L0Pt  =(AtV_1A)  1AtV'1  [1.17] 

If  further  we  may  assume  that  the  errors  on  the  slope  measurements  have  equal  variance  and  are 
uncoirelated  (i.e.  V  =  cr2I )  we  obtain  a  standard  least-squares  solution  - 
Lop^fA^A1  [1.18] 


We  have  implemented  modal  wavefront  estimation  in  our  model  using  both  the  standard 
least-squares  solution  (eq.  [1.17])  and  the  optimal  MAP  estimator  (eq.  [1.18]).  The  current  version 
of  the  software  uses  the  well-known  Zemike  circular  polynomials  [5]  as  a  basis  set. 

1.5  Adaptive  mirror 

The  simulated  wavefront  correction  device  was  based  on  an  independent  tip-tilt  mirror  and  a  17 
element  bimorph  mirror.  The  influence  functions  used  to  represent  the  action  of  the  mirror 
electrodes  follow  the  theoretical  model  of  Halevi  [6],  This  model  assumes  that  the  mirror  is 
composed  of  one  piezoelectric  plate  glued  with  a  plate  of  mirror  glass  activated  by  a  distribution 
of  electrodes  placed  on  the  reverse  side  of  the  mirror.  The  device  has  a  rectangular  shape  of  size 
ab  m2.  In  such  a  model,  the  vertical  displacement  of  the  mirror  W(x,y)  and  the  voltage 
distribution  V(x,y)  must  satisfy  the  following  differential  equation  - 
V*W+AVV2  =  0  [1.19] 


The  coefficient  A  is  a  constant  proportional  to  dutm  2,  where  du  is  a  coefficient  of  the 

piezoelectric  tensor1  and  tm  is  the  thickness  of  the  mirror.  Typically  A  is  of  magnitude  10'9  m‘]V" 
.  The  model  assumes  that  the  mirror  is  constrained  at  the  edges  (i.e.  W=0  at  the  edges).  Then, 
by  representing  the  voltage  distribution  as  a  double  Fourier  series,  and  inserting  it  into  eq. 
[1.19],  it  follows  that  the  vertical  displacement  can  be  expressed  [6]  as  - 


W(x,y)=ZZ 

n=  1  m=  1 


v-  - 


[1.20] 


As  this  model  of  the  mirror  was  square,  we  use  a  circular  mask  to  match  the  mirror  to  the 
conjugate  image  of  the  pupil  (therefore  constraining  it  at  the  edge  of  the  mask).  The  distribution 
of  the  electrodes  was  chosen  as  a  double  circular  ring.  The  model  assumes  that  the  mirror  is 
linear  in  it’s  response  and  that  it  exhibits  no  hysteresis,  assumptions  which  are  not  exactly 
fulfilled  in  practice  but  which  nonetheless  do  not  generally  exhibit  a  major  influence  on 
performance.  The  time  response  of  the  mirror  is  not  modelled  but  it’s  overall  effect  on  imaging 
can  effectively  be  included  by  introducing  an  overall  time-delay  between  calculation  of  the 
required  drive  voltages  and  their  application  to  the  actuators  and  this  is  included  in  the 
simulation. 

1.6  Closed-loop  control  algorithm 

Our  computer  model  allows  the  use  of  either  an  integral  control  algorithm  or  a  Smith 
compensator  [7],  In  both  cases,  the  control  system  applies  a  correction  to  the  bimorph  mirror 
device  every  x  seconds  and  the  mirror  profile  is  maintained  during  this  interval  r  until  the  next 
correction  is  made.  The  correction  to  the  mirror  comprises  an  adjustment  to  the  voltages  applied 
to  the  actuators  and  is  calculated  on  the  basis  of  the  difference  between  the  residual  wavefront  as 
measured  by  the  wavefront  sensor  and  the  profile  currently  existing  on  the  adaptive  mirror. 

Let  the  set  of  voltages  yk(n(r-l)  )  be  the  voltages  existing  on  the  actuators  at  time  n(z- 1) . 
The  voltages  applied  to  the  actuators  at  time  nz ,  yk(nr)  are  given  by  the  integral  recursion  - 
y(nz)  =  y((n-l)z)  +  Kczu(nz)  [1-21] 

where  the  set  of  modification  voltages  u(ni)  are  calculated  from  the  present  slope  data  provided 
by  the  wavefront  sensor  and  Kcis  the  gain  of  the  integral  controller.  The  time  delay  r  and  the 

gain  of  the  integral  controller  may  be  adjusted  by  the  user. 

The  Smith  controller  was  proposed  to  compensate  for  the  problem  of  time  delay  in  automatic 
control  systems  in  1956  by  O.J.  Smith.  At  that  time,  digital  electronics  was  not  sufficiently 
developed  to  allow  it's  easy  implementation.  It  is  now,  however,  a  well-used  method  in  control 
systems  as  a  device  primarily  added  to  the  integral  controller  to  increase  stability.  The  approach 
consists  in  adding  a  feed-back  loop,  which  simulates  the  behaviour  of  the  plant  but  without  the 
correction  time  delays  resulting  from  finite  calculation  and  response  time  during  the  data 
processing  and  control  actuation.  If  we  let  the  set  of  voltages  y(  n(T-l)  )  be  the  voltages  existing  on 
the  actuators  at  time  n(T-l),  for  the  Smith  compensator  the  voltages  applied  to  the  actuators  at  time 
nT,  y(nT)  are  given  by  the  integral  recursion  - 


y{nT)  = 


KT 

\-KT 


u(nT )  + 


yikn-\)T) 
1  -KT 


-lll—y((n-2)T) 
1  -KT 


[  1.22] 


where  the  set  of  modification  voltages  u(nT)  are  calculated  from  the  present  slope  data  provided 
by  the  wavefront  sensor  and  K  is  the  gain  of  the  controller.  The  time  delay  T  and  the  gain  of  the 
integral  controller  may  be  adjusted  by  the  user. 


The  Smith  controller  reduces  then  to  a  second  order  algorithm  as  opposed  to  the  first  order 
algorithm  of  the  integral  controller. 

1.7  Graphical  User  Interface  (GUI) 

The  aim  of  the  GUI  is  to  effect  an  easy  and  convenient  interface  for  the  user  to  carry  out 
simulations  of  adaptive  optics  imaging  scenarios.  A  related  aim  was  to  construct  the  interface  in 
such  a  way  as  to  reduce  the  possibility  of  unintentional  input/errors.  A  GUI  has  been  developed 
which  will  allow  the  input  of  key  simulation  parameters,  allow  the  easy  visualisation  and  storage  of 
results.  The  GUI  is  effective  in  its  basic  aims  though  there  is  scope  for  more  flexibility. 

1.8  Current  restrictions  on  the  adaptive  optics  simulation 

In  it’s  current  form,  there  are  certain  restrictions  on  the  type  of  adaptive  optics  system  and  the  range 
of  operating  conditions  that  can  be  modelled. 

•  The  program  considers  only  a  Shack-Hartmann  sensor  with  52  active  sub-apertures.  Note  that 
the  current  version  does,  however,  allow  variation  in  the  sub-aperture  size  expressed  in  terms  of 
the  Fried  parameter  r0 .  Critically,  this  does  allow  the  subaperture  size  to  r0  ratio  to  be  varied,  a 
parameter  which  is  known  to  significantly  affect  low  light  level  performance  [3,8], 

•  The  basis  matrix  for  estimation  of  atmospheric  wavefronts  is  modelled  in  terms  of  the  Zemike 
polynomials  only.  This  is  not  a  serious  restriction  but  a  future  version  should  allow  the  use  of 
the  Karhunen-Loeve  basis. 

•  The  adaptive  mirror  is  a  bimorph  device  with  17  actuators.  There  is  no  facility  to  consider  other 
types  of  correction  device. 

•  Although  a  variety  of  closed-loop  control  algorithms  are  possible,  in  principle,  the  current 
version  permits  an  integral  controller  or  a  Smith  compensator  only. 

•  Basic  program  output  is  the  time-averaged,  compensated  Strehl  ratio  in  the  image  plane  as  a 
function  of  the  controller  gain.  This  is  a  satisfactory  measure  of  the  basic  performance.  Program 
output  for  intermediate  stages  such  as  random  phase  screens  and  Hartmann  data  and  image 
plane  data  is  possible  but  can  not  currently  be  effected  from  the  GUI. 

1.9  Compiling  and  executing  the  simulation  code. 

The  simulation  source  code  consists  of  nearly  10  000  lines  of  C.  A  number  of  data  structures  are 
also  necessary  for  the  program  to  run.  If  at  least  one  pre-calculated  set  of  time-evolving  wavefronts 
are  included  this  occupies  approximately  100  Mbytes. 

•  The  source  code  files  (which  are  too  numerous  to  list)  can  be  transferred  on  request  to  users  as 
required  by  ftp. 

•  All  source  code  files  should  be  placed  in  a  single  directory. 

•  To  compile  the  program,  an  ANSI  C  compiler  is  required.  A  text  file prbc.txt  contains  the 
necessary  commands  to  compile  and  link  the  code. 

1.10  Using  the  proaram 

Within  the  constraints  specified  in  section  1 .8,  the  user  may  select  a  number  of  observational  and 
system  parameters  which  may  be  expected  to  influence  the  overall  performance  of  the  AO  system 
using  a  GUI.  The  start-up  GUI  is  shown  overleaf- 


DISPLAY  DATA 


COMMAND  LINE  INTERFACE 


The  program  may  also  be  run  using  a  simple  command  line  menu  which  asks  the  user  to  specify  the 
conditions  of  the  simulation.  A  typical  example  is  given  in  figure  1  and  an  explanation  of  the 
choices  made  by  the  user  follows. 


»Do  you  want  to  use  pre-set  values  for  the  telescope  diameter  and 
turbulence  conditions  (i.e.  D,  rO,  v..)  ?  (y/n) 

This  is  an  important  choice  because  the  calculation  of  a  sequence  of  statistically  correct,  time- 
evolving  atmospheric  wavefronts  takes  considerable  computational  resources  and  a  long  time  to 
generate  (many  hours).  If  you  have  a  pre-calculated  set  of  wavefronts  corresponding  to  a  given 
telescope  diameter  D  and  Fried  parameter  r0  which  you  wish  to  use,  you  should  answer  n.  In  the 

example  in  figure  1,  the  user  has  answered  n. 

»Enter  filename  containing  this  information. 

Enter  the  filename  containing  the  telescope  diameter  and  turbulence  parameters  for  which  the  time- 
evolving  wavefront  is  already  pre-calculated. 

»Do  you  want  to  use  a  set  of  system  parameters  already  existing  ?  (y/n) 

Here,  the  user  may  decide  whether  to  simulate  an  imaging  scenario  for  which  the  relevant 
parameters  are  prespecified  in  a  file  or  alternatively  specify  new  system  parameters  which  will  be 
written  to  file  themselves.  In  the  example  in  figure  1,  the  user  has  answered  n. 

»Enter  file  name  for  saving  new  simulation  parameters  i 

An  appropriate  file-name  is  entered. 

»Infinite  or  finite  light  level?  (1  for  oo :  2  for  finite) 

The  object  being  imaged  must  be  specified  as  being  very  bright  (i.e.  providing  an  essentially  x 
light  level  at  the  wavefront  sensor)  or  a  finite  light  level.  In  the  example  in  figure  2,  the  user  has 
specified  a  finite  light  level. 

»Enter  visual  magnitude  of  the  natural  guide  star. 

Give  the  visual  magnitude  of  the  natural  guide  star  being  used  for  wavefront  sensing 

»Enter  sensitive  bandwidth  of  wavefront  sensor  camera  AA  (nm) 

The  user  must  specify  the  optical  wavelength  region  in  nanometres  to  which  the  wavefront  sensor 
camera  is  sensitive. 

»Enter  CCD  quantum  efficiency  (0  and  1) 

Specify  the  quantum  efficiency  of  the  CCD  camera 

»Specify  level  of  CCD  noise  (rms  photoelectrons/per  pixel/frame) 

The  user  must  specify  the  noise  performance  of  the  CCD  camera  in  terms  of  the  rms  photoelectrons 
generated  per  pixel  per  camera  frame. 

»Enter  minimum  gain  for  integral  controller:  kmin _ 

The  output  of  the  AO  simulation  is  the  time-averaged,  compensated  Strehl  ratio  as  a  function  of  the 
gain  of  the  integral  controller.  Specify  the  minimum  gain  of  interest. 

»Enter  maximum  gain  for  integral  controller:  kmax _ 

The  output  of  the  AO  simulation  is  the  time-averaged,  compensated  Strehl  ratio  as  a  function  of  the 
gain  of  the  integral  controller.  Specify  the  maximum  gain  of  interest. 

»Enter  overall  time  delay  between  wavefront  measurement  and  correction  on  mirror  (ms) 

Any  real  adaptive  optics  system  operating  in  closed-loop  will  introduce  a  delay  between  calculating 
the  current  form  of  the  residual  wavefront  and  applying  an  appropriate  correction  to  the  mirror. 
Specify  this  delay  in  ms. 

»Enter  numeric  label  for  storing  result  of  simulation _ _ _ _ 

This  allows  the  user  to  label  the  simulation  results  file  via  a  parameter  c_offset  which  is  set 


»Enter  any  5  digit  random  number  (for  noise  simulation) _ 

This  ensures  that  each  simulation  will  begin  with  a  different  random  number  seed. 


»Enter  directory  label  _ _ _ _ _ 

This  allows  the  user  to  specify  the  directory  in  which  to  place  the  simulation  results  file.  The 
parameter  dirhold 

The  simulation  is  carried  out  and  the  results  are  stored  in  a  file  named  — /*/zern_result*.d.  The  first 
asterisk  corresponds  to  the  directory  label  dir_hold  and  the  second  to  the  numeric  label  c_offset.  In 
other  words,  if  the  program  is  executed  from  the  directory  “current_working_directory”  choosing 
dir  hold  =  myplace  and  c  offset  =  3  will  cause  the  output  file  to  be  written  as  - 
—/current  working  directorv/mvplace/zern  result3.d  .This  file  contains  the  compensated  Strehl 
ratios  achieved  for  each  value  of  the  integral  gain 


2  Low-light  level  simulations  -  some  preliminary  results 


In  this  section,  we  present  results  of  the  performance  of  a  Shack-Hartmann  based  adaptive 
optics  system  operating  under  low  photon  flux  conditions.  The  new  generation  of  large,  ground- 
based  telescopes  in  the  5-10  m  diameter  range  which  are  currently  being  built  and  commissioned  in 
many  parts  of  the  world  (e.g.  ALFA,  GEMINI,  SUBARU)  will  all  be  equipped  with  adaptive  optics 
systems  employing  laser-guide  stars.  We  thus  stress  that  the  relevance  of  the  presented  results  lies 
more  in  the  use  of  smaller  ground-based  telescopes  (possibly  several  small  telescopes  employed  in 
an  interferometric  configuration)  where  more  modest  D/rO  values  prevail  and  laser-guide  star 
technology  is  either  not  available  or  undesirable. 

The  series  of  simulations  assumed  the  following  fixed  conditions  - 

1 .  D/rO  =  5.  Telescope  diameter  lm  and  rO  =  20cm. 

2.  Average  wind  speed  of  20  m/s.  For  a  Fried  parameter  of  r0=20cm,  this  gives  a  Greenwood 
time-delay  of  3.15  ms. 

3.  A  Shack-Hartmann  wavefront  sensor  having  52  active  sub-apertures. 

4.  A  CCD  camera  of  fixed  quantum  efficiency. 

5.  A  17  actuator  bimorph  mirror. 

The  following  conditions  were  varied  - 
The  incident  photoelectron  flux. 

The  read-noise  in  photoelectrons  per  pixel  per  frame  on  the  CCD  camera. 

The  time-lag  r  between  wavefront  sensing/estimation  and  application  of  the  actuator  signals 
The  gain  on  the  closed-loop  control  algorithm  (either  the  integral  controller  or  the  Smith 
Compensator)  was  varied. 

The  aim  was  to  study  the  mean  compensated  Strehl  ratio  of  the  simulated  AO  system  as 
these  parameters  varied.  In  particular,  the  optimisation  of  the  closed-loop  control  gain  was  studied 
for  several  different  signal-noise/correction  time-delay  scenarios. 

Figure  1  shows  the  mean-compensated  Strehl  ratio  as  a  function  of  controller  gain.  In  this  study,  the 
incident  flux  was  high  (i.e.  assumed  to  be  infinite).  The  control  algorithm  was  the  integral 
controller  and  the  correction  time-lag  x  was  set  =1.5  ms.  This  time-lag  is  only  approximately  1  h 
of  the  Greenwood  time  delay  and  the  system  is  thus  fast  enough  to  accurately  track  the  atmospheric 
fluctuations.  This  rather  idealised  case  allowed  us  to  test  the  limiting  behaviour  of  the  integral 
controller.  We  note  an  optimum  gain  around  400  (this  corresponds  to  a  linearisation  factor 
K,r  ~  1).  As  we  go  to  higher  values  of  gain  the  system  becomes  unstable  due  to  the  spatial 
modelling  error  and  the  fact  that  the  time  delay,  although  small,  is  still  finite. 

Figure  2  shows  the  mean-compensated  Strehl  ratio  as  a  function  of  controller  gain.  In  this  study,  the 
incident  flux  was  low,  corresponding  to  a  mean  value  of  60  photoelectrons/sub-aperture.  The 
control  algorithm  was  again  the  integral  controller  and  the  correction  time-lag  r  =  1.5  ms.  A  read- 
noise  level  of  1.5  photoelectrons  per  pixel  per  frame  (state-of-the-art  devices  currently  achieve 


values  of  -  3  )  was  also  included.  The  main  features  of  figure  2  are  largely  in  agreement  with 
intuition.  Overall  compensated  Strehl  is  reduced  and  the  gain  curve  is  more  peaked  with  a 
maximum  occurring  around  KI  ~  3  5  0 .  The  rapid  decrease  in  performance  above  the  optimum 
reflects  the  tendency  of  the  integral  controller  to  amplify  the  noise. 

Figure  3  was  obtained  under  identical  conditions  to  figure  2  except  that  the  Smith  controller  was 
implemented.  In  this  case,  we  see  that  the  optimum  gain  is  somewhat  lower  and  overall 
performance  slightly  poorer  for  high  gain.  As  the  correction  time  delay  is  much  smaller  than  the 
Greenwood  time,  we  do  not  expect  any  significant  improvement  from  the  elimination  of  the  time 
delay  and  the  poor  noise  propagation  for  high  gain  results  in  a  slightly  lower  performance. 

Figures  4  and  5  show  the  mean-compensated  Strehl  ratio  as  a  function  of  controller  gain  for 
integral  and  Smith  controllers  with  a  large  correction  time  delay  of  8ms  (>  2.5  times  the 
Greenwood  time).  The  flux  was  100  photoelectrons  per  sub-aperture  and  the  read-noise  at  a  level  of 
1.5  photoelectrons/sub-aperture/ffame.  Despite  the  increase  in  signal-noise,  figures  4  and  5  show  a 
decrease  in  the  overall  compensation  as  compared  to  figures  2  and  3.  This  results  from  the 
significant  correction  time  delay.  We  note  that  in  this  case,  the  Smith  compensator  performs 
slightly  better  than  the  integral  controller  and  is  stable  across  a  larger  range  of  gains. 

These  studies  have  demonstrated  the  ability  of  the  program  to  simulate  AO  systems 
operating  in  low  light  level  regimes.  Significant  improvements  in  Strehl  ratio  over  the 
uncompensated  state  (up  to  S=  0. 13)  can  be  achieved  even  at  flux  levels  of  60  photoelectrons/sub- 
aperture  using  AO  compensation  which  uses  a  very  low-noise  CCD  camera  (1.5  electrons 
rms/pixel/frame)  in  the  wavefront  sensing  arm.  We  have  also  demonstrated  the  ability  of  the  Smith 
controller  to  partially  compensate  the  correction  time-delay  and  yield  improved  performance  over 
the  standard  integral  controller.  A  critical  aspect  in  these  low  flux  studies  was  the  read-noise 
performance  of  the  simulated  CCD  camera  -  a  factor  of  3  increase  in  the  read-noise  under  the  same 
conditions  was  sufficient  to  effectively  reduce  AO  compensation  to  a  performance  level  similar  to 
uncompensated  imaging.  So,  while  we  stress  that  greater  sensitivity  in  the  wavefront  sensing  can 
lead  to  enhanced  performance  and  can  be  achieved  via  several  strategies  that  have  been  studied  in 
related  work  (namely,  matched  filtering  techniques  for  centroiding  [4],  optimal  Bayesian  wavefront 
estimators  [9]  and  also  the  crucial  choice  of  the  Hartmann  sub-aperture  size  to  match  the  prevailing 
signal-noise  level  [3],  the  importance  of  using  low  noise  imaging  devices  in  the  wavefront  sensing 
arm  must  be  emphasised.  As  current  state-of-the  art  CCD  devices  yield  read-noises  of  the  order  of 
3  electrons  rms/pixel/frame,  photon  counting  devices  which  have  effectively  no  noise  (but  currently 
considerably  lower  quantum  efficiency)  may  prove  to  be  the  right  instrument  of  choice  for  low 
light  level  AO.  A  simulated  photon  counting  device  having  20  %  quantum  efficiency  but  no  read- 
noise  yields  a  broadly  similar  performance  to  that  given  in  figure  1.  A  systematic  comparison 
between  such  devices,  however,  remains  to  be  carried  out. 


FIGURE  1:  MEAN  COMPENSATED  STREHL  RATIO  /  GAIN 
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FIGURE  4:  MEAN  COMPENSATED  STREHL  RATIO  /  GAIN 
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FIGURE  5:  MEAN  COMPENSATED  STREHL  RATIO  /  GAIN 
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Section  2:  Summary 


•  The  software  is  able  to  simulate  the  imaging  capacity  of  an  AO  system  viewing  an  object 
confined  to  a  region  within  a  single  isoplanatic  patch.  The  predominant  interest  is  in  using  this 
system  for  low  light  level  studies.  No  laser  guide  star  system  is  assumed  and  no  science-object 
is  necessary  to  assess  performance.  The  AO  system  consists  of  a  17  element  bimorph  adaptive 
mirror  coupled  to  a  Hartmann  sensor  employing  52  active  sub-apertures  and  a  low-noise  CCD 
camera  operating  in  closed-loop  with  an  integral  or  Smith  controller.  The  incident  flux  from  the 
wavefront  sensor  object,  turbulence  conditions,  telescope  diameter,  type  and  gain  of  the 
controller,  correction  time  delay  and  CCD  camera  noise  performance  may  all  be  adjusted  by 
the  user.  Preliminary  results  suggest  that  the  software  performs  correctly  and  have  enabled 
basic  low  light  level  imaging  performance  to  be  assessed  as  a  function  of  the  variable 
parameters. 

Future  improvements  could  be  made  in  two  main  areas- 

•  Incorporating  greater  flexibility  into  the  GUI  to  control  a  larger  range  of  study/observational 
condition. 

•  Addition  of  software  modules  to  include  different  adaptive  mirrors  or  correction  devices  (e  g. 
LCs),  different  imaging  devices  such  as  photon-counting  cameras  and  more  sophisticated 
approaches  to  wavefront  estimation  and  control  (e.g.  Matched  filtering  techniques  for  centroid 
estimation  and  predictive  control  algorithms). 
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Abstract 

Common  wavefront  sensors  such  as  the  Hartmann  or  curvature  sensor  provide  mea¬ 
surements  of  the  local  gradient  or  Laplacian  of  the  wavefront.  Expression  of  wavefronts 
in  terms  of  a  set  of  orthogonal  basis  functions  thus  generally  leads  to  a  linear  wavefront 
estimation  in  which  modal  cross-coupling  occurs.  Auxiliary  vector  functions  may  be  de¬ 
rived  which  effectively  restore  the  orthogonality  of  the  problem  and  enable  the  modes  of 
a  wavefront  to  be  independently  and  directly  projected  from  slope  measurements.  Using- 
variational  methods,  we  derive  the  necessary  and  sufficient  condition  for  these  auxiliary 
vector  functions  to  have  minimum  error  norm.  For  the  specific  case  of  a  slope-based  sensor 
and  a  basis  set  comprising  the  Zernike  circular  polynomials,  these  functions  are  precisely 
the  Gavrielicles  functions. 


Ill  problems  of  wavefront  estimation,  it  is  common  to  expand  the  wavefront  in  terms 
of  a  set  of  orthogonal  functions  or  modes  - 

4>(i)  =  J^o.tPkix)  (1.1) 

k- 1 

which  obey  the  orthogonality  relation  - 

I  Pi(x)Pj{x)d2.r  =  6ij  (1.2) 

Jd 

and  D  denotes  the  domain  of  integration. 

The  most  commonly  used  set  of  functions  is  the  Zernike  circular  polynomials.  Speak¬ 
ing  generally,  orthogonal  bases  are  desirable  because  they  allow  the  modal  coefficients  to 
be  evaluated  by  simple  integration  of  a  product  of  two  functions  over  the  domain,  using 
the  orthogonality  relation  ecp  1.2.  However,  since  Hartmann  sensors  (and  shearing  inter¬ 
ferometers)  provide  estimates  of  the  gradient  of  the  phase  rather  than  the  phase  itself,  the 
appropriate  model  is  - 

(V*(xj))±i2ab(VPk(z.j))  .  (1-3) 

k  =  L 

where  the  coordinates  Xj  denote  the  position  of  the  slope  measurements  within  the 
pupil.  In  this  case,  the  orthogonality  of  the  basis  cannot  be  exploited.  Modal  cross¬ 
coupling  will  occur  and  the  coefficients  must  be  obtained  by  solving  an  inverse  prob- 
lem/overdeterminecl  system  of  linear  equations.  Many  workers  have  examined  this  problem 
in  a  search  for  optimal  solutions  [1-6]. 

The  evaluation  of  modes  by  direct  integration  could  be  restored  if  we  can  derive  a 
set  of  auxiliary  vector  functions  F  JxJ  which  are  orthogonal  to  the  gradients  of  the  basts 
functions .  To  establish  the  necessary  conditions  for  such  a  set  to  exist,  we  consider  that 
these  functions  obey  the  relation  - 

V-F  i(x)  =  Pi(.v)  (1.4) 


Using  the  orthogonality  relation  eq.  1.2  in  eq.  1.1  and  substituting  eq.  1.4.  the  modal 
coefficients  are  given  by  - 

(H  =  /  <p(-L)PiU)d~*  =  /  6(x)V-F i(x)d-.r 
Jd  Jd 


This  niav  be  written  as  - 


/  V-(F iU)o(.r))d-.v-  I  Vo(r).  FJrfd-x 

Jd  Jd 


o:i  = 


Applying  the  divergence  theorem,  we  obtain  - 


=  <p  6{x)Fi(x)  ■  dl  —  /  nablo,(b{x)  ■  ¥ i{x_)d“x 


Clearly  if  the  required  set  of  functions  Fz(r )  satisfies  the  relations  - 

V  •  F;(;r)  -  Pt(x) 

F i(x)  -cU  =  0  (1.6) 

where  Iq  ( x )  •  cU  denotes  the  normal  component  of  the  function  to  the  closed  contour  C  of 
the  integration  domain  D,  we  then  have  - 

cH  =  -  [  Vd>(x)  •  F i(x)d-x  (1.7) 

■ID 

In  other  words,  we  may  evaluate  the  modes  by  direct  integration  as  required.  It  has  been 
shown  that  there  are,  in  fact,  a  number  of  possible  sets  of  vector  polynomials  which  may¬ 
be  used  in  eq.  1.7  [7-9]  and  there  is  no  unique  solution  to  the  problem  described  by  eq. 


In  the  remainder  of  this  letter,  we  will  address  the  question  of  how  to  derive  the 
optimum  set  of  vector  functions  Ft(x)-  By  optimum,  we  mean  that  set  of  vector  functions 
F i(¥)  which  will  give  the  minimum  error  or  noise  propagation  in  our  estimate  of  the 
wavefront  given  by  eq.  1.7. 

The  estimate  of  the  kth  modal  coeffficient  is  given  by  - 


aic  =  ~  ju(z)  ■  F k  (£)d2x 
Jd 


where  p.  (  x_)  —  V  <p  (rj  4-  u  ( r ) ,  and  v  ( rj  is  the  additive  noise  vector  at  x_.  If  we  may 
assume  a  noise  process  having  zero  mean  and  covariance  cr^b(x  —  x')  (which  is  reasonable 
for  Hartmann-type  sensors),  it  is  straightforward  to  show  that  the  variance  in  radians2 
associated  with  our  estimate  of  the  kth  mode  is 


N(Fk)  =  <rl  /  |F*  (iiL)| -d-x 
Jd 


The  noise  propagator  is  thus  defined  as 


-V(Ffc)=  /  \Fk(£)\2d--r 


This  noise  propagator  thus  depends  only  on  the  particular  choice  of  the  vector  functions 
through  the  volume  defined  by  \F2k  (.r)|in  D.  The  ensemble- averaged  mean-square  error 
associated  with  the  estimator  is  then  - 

<'2  =  *;y>V(Fk>+  E  (4)  nil) 

i  =  2  k-M  +  L 

where  the  latter  term  simply  corresponds  to  unestimated  modes.  As  the  conditions  spec¬ 
ified  in  ecp  1.6  do  not  determine  a  unique  solution  for  F'k(x).  there  are  many  possible 
solutions  for  each  mode,  each  having  a  different  associated  noise  propagator. 

The  key  point  is  that  we  seek  to  minimise  the  noise  propagators  expressed  by  eq.  1.10 
subject  to  the  constraint  equation  and  boundary  condition  expressed  by  eq.  1.6.  This  is 
a  variational  problem.  Accordingly,  we  introduce  the  Lagrange  multiplier  function  \k(x) 
and  seek  to  minimise  an  objective  function  given  by  - 

Q  =  J  Fl(xyrz  +  j  \k(x)[V  ■  Ffc («)  -  Pk(x)]d2x  (1.12) 

We  now  take  first  variations  in  Q  with  respect  to  F&(r)  and  set  these  to  zero  to  find 
the  stationary  points  of  Q  - 

6Q  =  2  J  $Fk(x)  ■  Fjc(x_)cl"x  +  J  A fc(x)V  •  SF ic(x_)d" x  =  0  (1.13) 

Since  V  •  [Xk{x)Fk(x)}  =  A*(*)V  •  F*(z)  +  Fk(x)  •  VA*(x)  this  may  be  expressed  as  - 
SQ  =  2  j  SFk(x) •  F jfe (x)d~x -t-  j  V  ■[\k{x)5Fk{x)}c!2x- J  6Fk{x)- VXk{x)d2x  =0  (1.14) 

Applying  the  divergence  theorem  to  the  second  integral  in  ecp  1.14  and  rearranging  we 

obtain  -  ,. 

i[2Fk(x)  -  ■  SFk(x)d2x  +  j  Xk(x)SFk{x)  ■  dl_  =0  (1.15) 

The  boundary  condition  f  Fk(x)  ■  cl(_l).  =  0  from  eq.  1.6  ensures  that  the  last  integral  in 
eq.  1.15  is  identically  zero  and  we  arrive  at  the  simple  condition  for  the  Fj,.(r)  to  have 
minimum  error  norm  - 

FtU)  =  VA*(i)  (1.16) 

In  other  words,  the  Fjfe(:c)  are  an  irrotational  set  of  functions  -  i.e.  must  be  expressible 
as  the  gradient  of  an  associated  scalar  function.  Note  we  have  dropped  the  factor  of  2  as 
this  will  have  no  effect  on  the  functional  form  for  VA k{-L))-  Substituting  eq.  1.16  into  the 
constraint  eq  1.6  we  find  that  the  required  scalar  field  A*  (x)  is  determined  by  the  Poisson 
equation  with  Neumann  boundary  condition  - 

V-Ajt  (x)  =  Pk  (x)  x^D 


V  At  (r)  ■  n(r)  =  0  r  t  C 


(1.17) 


This  is  the  main  result  of  our  analysis.  The  Poisson  equation  with  A eumann  boundary 
conditions  has  a  unique  solution  and  the  minim  um  error  norm,  modal  projector  functions 
can  be  obtained  by  solving  eq.  1.17  for  the  given  set  of  basis  functions.  Pkia)*  Such  a 
set.  of  functions  was.  in  fact,  explicitly  derived  in  analytic  form  for  the  Zernike  circular 
polynomial  basis  by  Gavrielides  in  1982  [10]  although  the  noise  propagation  behaviour  of 
these  functions  was  not  addressed  in  his  original  work.  For  completeness,  their  analytic 
form  is  given  here  - 

The  components  of  the  Gavrielides  vector  polynomials  are  defined  in  polar  coordinates 
as  follows  - 

For  ni  ~  0 

F'r  ~  +  1  )T™  (r)co$m8  if  i  even 

7T 

F*  =  —  \/2(n  -f  1  )T™{r).$inm9  if  i  odd 


Fjj  =  1  >/2(n  +  l)Q™(r)[— m-sinmO]  if  i  even 
Fi  =  1  v/2(n  +  1)Q™ (r)[mcosm6]  if  i  odd 

7T 

where 


n  —  m  -  i  1 

1  ^  Cr(,)(n-2,  +  2)[r-^-r-^1] 

nK']  4  2-i  (ii±^  -  .s  +  1)  s  +  1) 


5=0 


« (r)  ’  4m  I! 


n  —  m 

1  ^  C?(s)  [(n  —  2s  +  2)  rn;— 1  -  mrn-t+l _ 


5=0 


(S^a_.s  +  1)  (a^a_3  +  i) 


and  where  C‘”l(s) 


-!*('»- 


•s) 


For  m=0 


1  , - -v^  (-l)'(n  - 


F:.  =  ^vrTTZ 


sir 


fz'od  [($-*)!] 3  ("-2* +  2) 


4  =  0  .  (1.18) 

Given  the  inherent  simphcity  of  modal  projection  from  slope  measurements  according 
to  eq.  1.7,  it  is  somewhat,  surprising  that  these  functions  have  not  seen  more  widespread  use 


in  wavefront  sensing  and  estimation  problems.  We  anticipate  that  the  knowledge  of  their 
minimum  noise  propagation  behaviour  will  make  them  a  simple  and  attractive  approach 
to  wavefront  estimation  from  slope  measurements  and  lead  to  their  increasing  use. 

If,  is  interesting  to  note  that  the  problem  we  have  discussed  of  obtaining  a  set  of 
auxiliary  vector  functions  which  gives  minimum  wavefront  error  norm  has  a  close  analogy 
with  minimum  physical  principles.  Eq.1.17  may  in  fact  be  identically  used  to  describe  the 
equilibrium  shape  of  a  membrane  (with  fixed  boundary  described  by  the  closed  contour 
C)  over  which  a  force  P&  (r)  acts  to  produce  vertical  displacement  X(x_ ).  In  this  case, 
f  F'l  (x)  cl2x  represents  the  potential  energy  stored  in  the  membrane  and  will  be  equal  to 
the  work  done  in  bringing  the  membrane  to  it  s  equilibrium  position  (/P&  (#)  X(x)d~x_  ) 
only  in  the  case  of  conservative  fields  -  i.e.  when  Fk,  {x)  =  VA*(:r)  [11]- 
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